Natural vectors of Plasmodium knowlesi and other primate, avian and ungulate malaria parasites in Narathiwat Province, Southern Thailand

To date, four species of simian malaria parasites including Plasmodium knowlesi, P. cynomolgi, P. inui and P. fieldi have been incriminated in human infections in Thailand. Although the prevalence of malaria in macaque natural hosts has been investigated, their vectors remain unknown in this country. Herein, we performed a survey of Anopheles mosquitoes during rainy and dry seasons in Narathiwat Province, Southern Thailand. Altogether 367 Anopheles mosquitoes were captured for 40 nights during 18:00 to 06:00 h by using human-landing catches. Based on morphological and molecular identification, species composition comprised An. maculatus (37.06%), An. barbirostris s.l. (31.34%), An. latens (17.71%), An. introlatus (10.08%) and others (3.81%) including An. umbrosus s.l., An. minimus, An. hyrcanus s.l., An. aconitus, An. macarthuri and An. kochi. Analyses of individual mosquitoes by PCR, sequencing and phylogenetic inference of the mitochondrial cytochrome genes of both malaria parasites and mosquitoes have revealed that the salivary gland samples of An. latens harbored P. knowlesi (n = 1), P. inui (n = 2), P. fieldi (n = 1), P. coatneyi (n = 1), P. hylobati (n = 1) and an unnamed Plasmodium species known to infect both long-tailed and pig-tailed macaques (n = 2). The salivary glands of An. introlatus possessed P. cynomolgi (n = 1), P. inui (n = 1), P. hylobati (n = 1) and coexistence of P. knowlesi and P. inui (n = 1). An avian malaria parasite P. juxtanucleare has been identified in the salivary gland sample of An. latens. Three other distinct lineages of Plasmodium with phylogenetic affinity to avian malaria species were detected in An. latens, An. introlatus and An. macarthuri. Interestingly, the salivary gland sample of An. maculatus contained P. caprae, an ungulate malaria parasite known to infect domestic goats. Most infected mosquitoes harbored multiclonal Plasmodium infections. All Plasmodium-infected mosquitoes were captured during the first quarter of the night and predominantly occurred during rainy season. Since simian malaria in humans has a wide geographic distribution in Thailand, further studies in other endemic areas of the country are mandatory for understanding transmission and prevention of zoonotic malaria.


Results
Distribution and abundance of Anopheles species. Altogether, 367 female Anopheles mosquitoes were caught during rainy and dry seasons in 2018 and 2019. The total number of Anopheles collected from Sukhirin and Waeng Districts in Narathiwat Province during the raining seasons was 2.7 times more than those captured during the dry seasons ( Fig. 1, Table 1). Based on taxonomic keys and molecular analysis of Anopheles 33 , the predominant species belonged to An. maculatus, An. barbirostris s.l., An. latens and An. introlatus accounting for 96.19% of all Anopheles mosquitoes while those identified as An. umbrosus s.l., An. hyrcanus s.l., An. minimus, An. aconitus and An. kochi were sporadically found, ranging from 1 to 4 specimens for each species (Table 1). The biting rates for An. maculatus, An. barbirostris s.l., An. latens and An. introlatus were 3.400, 2.875, 1.625 and 0.925 mosquito collected/night/collector, respectively. The number of Anopheles mosquitoes captured for each corresponding period of the night between 2018 and 2019 seemed to show a similar trend. The overall number of mosquito collection peaked at 20:00 and 21:00 h whereas remarkably fewer numbers of mosquitoes could be caught after midnight, especially after 2:00 h until daybreak ( Fig. 2A). In the mosquito collection sites, the level of temperature peaked during 21:00 and 23:00 h whereas humidity gradually rose from dusk till dawn. The abundance of mosquitoes during the first quarter of the night (18:00 and 21:00 h) seemed to be positively correlated with the environmental temperature (Pearson r = 0.999, p = 0.03) and roughly correlated with the levels of humidity (Pearson r = 0.996, p = 0.05) ( Fig. 2A). The number of An. latens and An. introlatus seemed to be early feeders from 18:00 and 19:00 h whereas peak feeding time of An. barbirostris s.l. was between 20:00 and 21:00 h. Meanwhile, An. maculatus was most abundant from 21:00 to 23:00 h. Despite the decline in the number of mosquitoes after midnight, both An. maculatus and An. barbirostris s.l. could be sparsely caught before dawn (Fig. 2B).
Identification of malaria parasites. Of 367 Anopheles mosquitoes analyzed, primary PCR targeting the mitochondrial cytochrome c oxidase subunit 1 (cox1) gene of Plasmodium revealed positive results in 19 samples, accounting for 5.2% of total samples. Most Plasmodium-positive salivary glands (14 of 19 samples, 73.7%) were obtained during the rainy seasons. Species-specific PCR targeting human and simian Plasmodium species could detect 11 mosquitoes bearing simian malarial DNA in their salivary glands including P. inui (n = 5), P. fieldi (n = 2), P. knowlesi (n = 1), P. coatneyi (n = 1), P. cynomolgi (n = 1) and co-existence of P. knowlesi and P. inui (n = 1). DNA sequencing of recombinant plasmid clones from Plasmodium-positive specimens has reaffirmed the species of almost all simian malaria parasites except two P. inui-positive samples (HBT177 and HBT353) whose sequences www.nature.com/scientificreports/ virtually belonged to P. hylobati ( Table 2). The species of Plasmodium in the remaining eight mosquitoes could not be determined by species-specific nested PCR. Sequences from recombinant plasmid clones of these unassigned samples displayed six different Plasmodium species/lineages by phylogenetic analysis. Of these, two mosquitoes (HBT181 and HBT368) harbored Plasmodium species that had phylogenetic affinity to an unnamed simian malaria parasite known to infect long-tailed and pig-tailed macaques in Malaysian Borneo (GenBank accession no. KJ569860) (Fig. 3) 34 . One sample yielded the cox1 sequence belonging to P. juxtanucleare with 99.79% sequence identity. Four other mosquitoes (HBT253, HBT329, HBT340 and HBT341), each containing 2 to 3 distinct cox1 alleles, seemed to diverge from P. circumflexum and were placed into three distinct phylogenetic clades. The p-distance (d ± S.E.) between the cox1 and its flanking sequences of P. circumflexum and those of the

Diversity of simian Plasmodium species from mosquitoes, humans and macaques in Thailand.
To determine allelic variation in the cox1 and its flanking sequences of simian malaria parasites among different hosts, the sequences obtained from the salivary glands of mosquitoes in this study were compared with the corresponding gene sequences previously reported from macaques and humans in Thailand 12,13 . Results revealed that the cox1 and its flanking sequences of P. knowlesi, P. cynomolgi and P. fieldi were different among

Molecular identification of Anopheles species.
The species of all Plasmodium-positive mosquitoes were determined by sequencing of the PCR-amplified 2140 bp region of Anopheles mitochondrial genes encompassing cox1 and cox2. Due to the lack of representative cox2 sequences of some members in the Leucosphyrus Group, species assignment was determined mainly from the cox1 phylogeny (Fig. 5A). Results revealed that An. introlatus carried 4 species of primate malaria parasites including P. knowlesi, P. cynomolgi, P. inui and P. hylobati, and a plausible novel species (HBT341) of an avian malaria parasite related to P. circumflexum.
In the salivary glands of An. latens, six species of primate malaria parasites were identified including P. knowlesi, P. inui, P. fieldi, P. hylobati, P. coatneyi and an unnamed macaque malaria parasite Plasmodium sp. (GenBank accession no. KJ569860). Furthermore, P. juxtanucleare and a plausible novel species of avian plasmodia were also found in An. latens (Table 2). It is noteworthy that P. caprae was detected in the salivary glands of An. maculatus in this study. Although the sequence spanning cox1 and cox2 of mosquito HBT329 bearing an unknown Plasmodium sp. related to P. circumflexum could not be assigned from BOLD database 36 , it has been identified as An. macarthuri from phylogenetic analysis inferred from available reference cox1 sequences of Leucosphyrus mosquitoes spanning the 250-bp fragments (Fig. 6). Likewise, the other unassigned species of mosquito (HBT165) without Plasmodium infection also belonged to An. macarthuri based on phylogenetic analysis (Fig. 6). Meanwhile, the tree constructed from the cox2 locus per se displayed concordant phylogenetic affinity of mosquitoes with that inferred from the cox1 locus (Fig. 5B).

Discussion
In Thailand, the past couple of decades have seen a dramatic decline in the number of falciparum malaria patients while a relative increase in the proportion of vivax malaria has been envisaged, suggesting that currently applied control measures seem to be less effective against non-falciparum infections [10][11][12][13]37,38 . Likewise, an increased prevalence of patients infected with P. knowlesi has been observed during the past decade, especially among those who resided in areas where infected domesticated or wild macaques co-existed [10][11][12][13]37 . In 2022, a total of    40 . However, differentiation between P. inui and P. hylobati could be resolved and correctly assigned by analysis of the cox1 sequences as shown in our previous study 13 . It is noteworthy that P. hylobati has not been detected among macaques in Thailand 13,31,32 while the prevalence of this simian malaria species in gibbons known as its natural hosts awaits further investigation. The presence of P. knowlesi DNA in the salivary glands of An. latens and An. introlatus in Narathiwat Province has implied that both mosquito members in the Leucosphyrus complex could serve as potential vectors of this simian malaria species in this province where both infected humans and macaques have been previously identified [10][11][12][13]31,32 . Mosquitoes belonging to the Leucosphyrus Subgroup comprise at least 13 species while  www.nature.com/scientificreports/ seven of these have been detected in Thailand 33 . The breeding places of these mosquitoes included ground pools, animal footprints, wheel-tracks and small shallow running streams in the forested ecotype with partial or heavilyshaded areas 33 [19][20][21] . Furthermore, An. latens could serve as a vector of P. coatneyi while the presence of P. coatneyi and P. hylobati DNA in the salivary glands of An. introlatus has suggested that they could transmit these simian malarias in the survey areas. In this study, clonal diversity in the cox1 sequences of most primate malaria parasites from infected mosquito salivary glands (Figs. 3, 4) and co-existence of different species of Plasmodium in a single mosquito (An. introlatus sample HBT206, Table 1) could probably stem from multiple feedings of the vectors on different monkeys infected with variant strains or different species of Plasmodium. Alternatively, multiple clones/species of malaria parasites in individual macaques could be prevalent, lending plausibility for the mosquitoes to acquire multiple clones/species of plasmodia upon single feedings. Our previous study has shown considerable clonal diversity of P. knowlesi and other primate malaria parasites in both long-tailed and pig-tailed macaques in Southern Thailand 32 ; thus, the latter scenario could not be excluded. Within mosquitoes, inter-species and intra-species competition could further influence subsequent transmission of malaria parasites to the vertebrate hosts 44,45 . Meanwhile, multiclonal and mixed species infections with simian malaria in humans occurred less frequently than in mosquitoes and macaques in Thailand 12,13,31,32 . Intriguingly, differential capability of adaptation to humans might occur among strains/lineages of primate Plasmodium species while genetic/phenotypic diversity of these parasites could contribute to a range of disease severity as observed in knowlesi malaria in humans albeit multiple host factors may also play crucial roles 4,5,8,[10][11][12][13]46,47 .
It has been shown that fluctuation in the levels of asexual parasitemia and gametocytemia occurred in experimental P. knowlesi infections in long-tailed macaques with peak parasitemia between 18:00 and 24:00 h and between 24:00 and 06:00 h, respectively 48 . The biting period of An. latens in Sarawak during 19:00 and 06:00 h seemed to coincide with the parasitemia of both asexual and sexual stages in peripheral blood of infected hosts; thus, enhancing transmission of simian malaria 16,19,23 . In field studies, An. latens and An. introlatus had peak biting activity between 19:00 and 20:00 h and 20:00 and 21:00 h, respectively 16,20,43 . Consistently, in this study simian Plasmodium-infected An. latens and An. introlatus were captured between 18:00 and 21:00 h ( Table 2), implying that exposure to mosquito bites in areas where infected macaques coexist in Southern Thailand during the first quarter of the night may pose a high risk for getting zoonotic malaria.
Phylogenetic trees inferred from the Plasmodium cox1 sequences derived from humans, macaques and Anopheles mosquitoes in Thailand have revealed genetic diversity across these hosts (Fig. 4). Despite limited number of available sequences, it is noteworthy that human isolates of P. knowlesi, P. cynomolgi, P. inui and P. fieldi differed from those derived from macaques and Anopheles mosquitoes, implying that the extent of genetic variation in natural populations of these primate malaria parasites could be more extensive. However, it seemed that a certain strain of P. inui derived from An. latens and An. introlatus in this study shared identical cox1 and its flanking sequences across 1504 bp with 11 previously characterized isolates from pig-tailed macaques from Choairong, Waeng and Sukhirin Districts in Narathiwat Province obtained in 2008, 2013, 2018 and 2020 13,32 . These isolates also possessed the same sequence as that of strain IC (leucosphyrus) derived from An. leucosphyrus caught from Negri Semilan in Malaysia in 1964 49 . Although identical mitochondrial cox1 sequences may not represent the same entire nuclear genomes, there remains a possibility that certain strains of P. inui could predominate and persist in natural transmission cycle that could probably stem from their adaptive capabilities to survive in diverse species of mosquito vectors and reservoir hosts.
Our sequence analysis of recombinant clones from 8 infected mosquitoes whose species of malaria parasites could not be assigned by PCR has revealed that 2 samples (HBT181 and HBT368) belonged to an unnamed Plasmodium known to infect both pig-tailed and long-tailed macaques in Sabah, northern Borneo (GenBank accession no. KJ569860) 34 . The placement of these isolates in the phylogenetic tree seemed to be closely related to P. coatneyi while the magnitude of evolutionary distance between the cox1 and its flanking sequences of these malaria species and that of P. coatneyi was comparable or greater than those between other primate Plasmodium species. Therefore, both isolates from An. latens (HBT181 and HBT368) and the previously described isolates (KJ569860) from Sabah 34 could belong to a distinct species. It is likely that more novel species of primate malaria parasites remain to be discovered. Meanwhile, Plasmodium species from An. maculatus (HBT314) displayed 99.77% to 99.92% sequence identity with that of P. caprae (GenBank accession no. LC090215). This ungulate malaria parasite has caused infections in domestic goats in Western Thailand and other countries 50  www.nature.com/scientificreports/ An. aconitus that have been recently reported to be implicated in transmission of P. caprae among goat farms in the western provinces of the country 51 . Apart from primate malaria, our analysis has identified An. latens as a potential vector of P. juxtanucleare, a predominant avian malaria parasite that infects chicken (Gallus gallus domesticus) and fighting cocks (Burmese red junglefowl, Gallus gallus spadiceus) in Thailand. Meanwhile, this avian malaria parasite has been known to infect other members of the Phasianidae 52-54 including short-crested flycatcher in Brazil 55 , turkey, parrot and chugar in Pakistan 56 , eared-pheasant in Japan 57 and black-footed penguins in South Africa 58 . Besides cosmopolitan in distribution, P. juxtanucleare has been detected among wild passerines in Brazil, suggesting spillover from poultry to free-living birds 59 . Although P. juxtanucleare is known to be mainly vectored by mosquitoes in the genus Culex (e.g. Cx. vishnui and Cx quinquefasciatus in The Philippines and Cx. saltanensis in Brazil) 60,61 , identification of An. latens as an additional vector in our study has expanded a range of mosquitoes capable of transmitting this avian malaria species. Herein, we also identified 3 Plasmodium lineages from An. latens and An. introlatus whose cox1 sequences displayed phylogenetic clustering within avian Plasmodium clades. Furthermore, the evolutionary distance (p-distance) between these lineages and a closely related P. circumflexum was comparable to or greater than those between known primate and avian malarial species, suggesting that these newly identified lineages could represent distinct avian Plasmodium species. Taken together, it is plausible that a wide range of mosquito vectors and diverse avian host species could enhance cosmopolitan distribution of avian malaria parasites.
Despite a wealth of knowledge about morphological differences among Anopheles vectors of malaria, molecular identification remains to be a useful method for fine resolution of mosquito species, especially those with minor structural differences. Although our sequence analysis encompassing near complete cox1 and partial cox2 sequences, An. latens and An. introlatus could be identified based on available cox1 reference sequences in the public database. The concordant topologies of phylogenetic trees inferred from either cox1 or cox2 sequences have provided an additional marker for speciation of these mosquito species (Fig. 5).
The midguts and ovaries of the mosquitoes have not been examined in this study. Although the presence of malarial DNA in the midguts may not replicate the results obtained from the mosquito's salivary gland samples, it can support the integrity of analysis from salivary gland specimens. Likewise, more information on the bionomics of the mosquitoes from examination of the ovary samples has not been available in this study. It is noteworthy that the presence of malarial DNA or sporozoites in the mosquito's salivary glands does not directly indicate its true vectorial status 1 . However, consistent identification of specific malarial DNA from the same species of mosquitoes, especially those from different geographic areas, seems to indicate that they could vector a given Plasmodium species 16 .
To the best of our knowledge, this is the first report on An. latens and An. introlatus as vectors of P. knowlesi in Thailand while both species could potentially transmit other simian and some avian plasmodia. Furthermore, An. maculatus has been identified to vector P. caprae that may parasitize domestic goats in the southern part of the country. Unraveling vectors of malaria and their bionomics could contribute to knowledge about disease transmission in malaria endemic areas that are important for prevention and control policy.  (Fig. 1). These locations were selected because (i) indigenous villagers residing in these areas were diagnosed to be infected with P. knowlesi [10][11][12][13] , (ii) both wild and domesticated pig-tailed and longtailed macaques were prevalent and known to be infected with several simian Plasmodium species 32 , (iii) both districts were situated along forest fringe with natural streams and rivers running through the areas that could be natural breeding places for Anopheles mosquitoes and (iv) rubber plantations were located in the districts where farmers and workers were frequently exposed to mosquito bites during the harvesting process. The collection sites had a tropical monsoon climate and were located near Hala-Bala tropical rain forest, covering an area of approximately 1.3 square kilometers across Yala and Narathiwat Provinces along Thailand-Malaysia border approximately between 5° 37′-6° 14′ N and 101° 8′ E-101° 51′ E. The elevation of Hala-Bala Forest ranged from 50 m to 1.5 km above mean sea level with annual mean temperature of 27.6 °C, annual average rainfall of 2560 mm, and relative humidity between 77 and 80%. The heavy rainy months of the year were from November to December and the relatively dry season from February to April. The number of days with rainfall (≥ 1.0 mm) occurred approximately 200 days per annum.

Collection of mosquitoes.
Mosquito collections were carried out every other night for 10 days during each of the four surveys in March and August 2018 and March and November 2019, covering relatively dry and rainy seasons of the study sites. Human-landing catches of mosquitoes were done outdoors. Each collection site comprised two teams, each consisting of two persons: one volunteer served as bait and the other as collector. The first team captured the mosquitoes from 18:00 to 24.00 h and the other from midnight until dawn (06:00 h). The collector captured the mosquitoes upon landing on the human bait by using a plastic collection tube. All captured mosquitoes were kept individually in a 1.5 ml sterile microtube with cap, placed on ice and transferred to 4 °C refrigerator for subsequent morphological identification. All volunteers and field staff were monitored for malaria infection during the field work and thereafter on a weekly basis for 2 months. Data analysis. Nucleotide sequences were aligned by using the codon-based option in MUSCLE program with manual adjustment 63 . For cox1 sequences of malaria parasites, known species or lineage of Plasmodium were used for comparative analysis which included the following GenBank accession nos.: P. falciparum (AJ276845), P. vivax (NC007243), P. malariae (AB354570), P. ovale curtisi (AB354571), P. ovale wallikeri (HQ712053), P. knowlesi (AB444106 and AY598141), P. cynomolgi (AB434919 and AB444126), P. inui (AB444110, AB444111, AB444114, AB444116 and AB444120), P. fieldi (AB354574 and KJ5698631), P. coatneyi (AB354575), P. fragile (AB444135 and AY722799), P. hylobati (AB354573), P. gonderi (AB434918), P. simium (AY800110), P. simiovale (AB434920), P. circumflexum (KY653762), P. juxtanucleare (AB250415), P. relictum (AY733090), P. elongatum (KY653802), P. caprae (LC090215) and P. sp. (KJ569860). The cox1 sequences of human-and macaque-derived P. knowlesi, P. cynomolgi, P. inui and P. fieldi in Thailand from our previous studies were included for comparison (GenBank accession nos. OL437412-OL437468 and OQ436033-OQ436048) 12,13 . For Anopheles cox1 and cox2 references, both complete and partial sequences of known species were deployed in the phylogenetic analysis: An. macarthuri (DQ897968 and DQ897969). All sites at which the alignment postulated a gap were eliminated in pairwise comparisons of the analysis. For evolutionary distance between sequences, p-distance was determined by dividing the number of nucleotide differences between two sequences by the total number of nucleotides compared and its standard error (d ± S.E.) between sequences was computed by 1000 bootstrap pseudoreplicates. Phylogenetic trees were constructed by using the maximum likelihood method with the best substitution model for the sequence data that yielded the minimum Bayesian Information Criterion (BIC) scores 64  www.nature.com/scientificreports/